Rasters in R

The first time you use a new version of R you will have to install all your packages. Try installing the raster and rgdal packages if you don’t have them already:

## Loading required package: raster

## Loading required package: sp

## Loading required package: rgdal

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /usr/share/proj
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

## Loading required package: lwgeom

## Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.8.0, PROJ 6.3.1

## Loading required package: sf

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1; sf_use_s2() is TRUE

## Loading required package: dplyr

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## Loading required package: osmdata

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

Now it is installed in your computer but isn’t loaded You will have to ‘import’ your library package every time you start up R.

library(raster)
library(rgdal)

Raster data is stored in a format similar to that of a matrix, except you have added data about georeferencing.

Let’s read in one raster and take a look

Examine Raster Data

notice that the working directory doesn’t need to be a path directly below your files of interest:

CWD_apr = raster('./data/raster/CWD/cwd2010apr.tif')

Get info about the raster object

CWD_apr

## class      : RasterLayer 
## dimensions : 1120, 872, 976640  (nrow, ncol, ncell)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : cwd2010apr.tif 
## names      : cwd2010apr 
## values     : 0, 150.9694  (min, max)

Plot the image:

plot(CWD_apr)

there are a variety of different plot options. Here are a few important ones

plot(CWD_apr,main = 'Title')

plot(CWD_apr,col= heat.colors(10), main='Water Deficit' ) # 10 designates the number of hex color codes that get generated

Now use ?plot to change col to a different preset color scheme

?plot

## Help on topic 'plot' was found in the following packages:
## 
##   Package               Library
##   terra                 /cloud/lib/x86_64-pc-linux-gnu-library/4.2
##   raster                /cloud/lib/x86_64-pc-linux-gnu-library/4.2
##   graphics              /opt/R/4.2.2/lib/R/library
##   base                  /opt/R/4.2.2/lib/R/library
##   sf                    /cloud/lib/x86_64-pc-linux-gnu-library/4.2
## 
## 
## Using the first match ...

Raster Math

Rasters like matrices can make basic math functions very easy. Note that rasters can be used in basic operations like CWD_apr *2 or CWD_apr^2.

However teo rasters with different extents, cell sizes, or projections CAN’T be used. You would need to reprojected and aligned first.

CWD_apr_sq = CWD_apr^2
plot(CWD_apr_sq, main='CWD Squared')

Now lets add two different rasters. Use the raster() function to read in the CWD for May

#CWD_may = raster('fill/in/the/path/file.tif')

CWD_A_M = CWD_apr + CWD_may
plot(CWD_A_M)

Raster Values

We can also directly set the cell values of rasters. For instance, it is often useful to have a raster filled with a fixed value, such as 0. In this example we will use the setValues() function to update the values of CWD_may to all be 0, and store it in a variable called zeros.

zeros = setValues(CWD_may,0)
plot(zeros)

Note that zeros has the same extent, resolution, and projection as CWD_may

print(zeros)

## class      : RasterLayer 
## dimensions : 1120, 872, 976640  (nrow, ncol, ncell)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : memory
## names      : cwd2010may 
## values     : 0, 0  (min, max)

print(CWD_may)

## class      : RasterLayer 
## dimensions : 1120, 872, 976640  (nrow, ncol, ncell)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## source     : cwd2010may.tif 
## names      : cwd2010may

We can also use boolean queries to update certain values. Here we update certain values of CWD_may to 5000 where CWD_may>150 is TRUE with the following syntax:

CWD_may[CWD_may>150] = 5000
plot(CWD_may)

Note that we have now changed the values of CWD_may… if we want to original values again, we will need to read them back in again with CWD_may = raster('./data/raster/CWD/cwd2010may.tif').

Challenge Question

In this section we will try to get the average water deficit for California for the whole year of 2010. To do this we will need to isolate the files we are interested in using - those ending with .tif.

First lets get a list of all files in our raster folder:

list.files(path="./data/raster/CWD/")

##  [1] "cwd2010apr.tfw"         "cwd2010apr.tif"         "cwd2010apr.tif.aux.xml"
##  [4] "cwd2010apr.tif.xml"     "cwd2010aug.tfw"         "cwd2010aug.tif"        
##  [7] "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml"     "cwd2010dec.tfw"        
## [10] "cwd2010dec.tif"         "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"    
## [13] "cwd2010feb.tfw"         "cwd2010feb.tif"         "cwd2010feb.tif.aux.xml"
## [16] "cwd2010feb.tif.xml"     "cwd2010jan.tfw"         "cwd2010jan.tif"        
## [19] "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml"     "cwd2010jul.tfw"        
## [22] "cwd2010jul.tif"         "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"    
## [25] "cwd2010jun.tfw"         "cwd2010jun.tif"         "cwd2010jun.tif.aux.xml"
## [28] "cwd2010jun.tif.xml"     "cwd2010mar.tfw"         "cwd2010mar.tif"        
## [31] "cwd2010may.tif"         "cwd2010nov.tif"         "cwd2010oct.tif"        
## [34] "cwd2010sep.tif"

Notice that .png and twf xml files are also listed. We can narrow this down by using a grep style pattern search.

Here we will only list files that have .tif in the name:

list.files(path="./data/raster/CWD/",pattern = '.tif')

##  [1] "cwd2010apr.tif"         "cwd2010apr.tif.aux.xml" "cwd2010apr.tif.xml"    
##  [4] "cwd2010aug.tif"         "cwd2010aug.tif.aux.xml" "cwd2010aug.tif.xml"    
##  [7] "cwd2010dec.tif"         "cwd2010dec.tif.aux.xml" "cwd2010dec.tif.xml"    
## [10] "cwd2010feb.tif"         "cwd2010feb.tif.aux.xml" "cwd2010feb.tif.xml"    
## [13] "cwd2010jan.tif"         "cwd2010jan.tif.aux.xml" "cwd2010jan.tif.xml"    
## [16] "cwd2010jul.tif"         "cwd2010jul.tif.aux.xml" "cwd2010jul.tif.xml"    
## [19] "cwd2010jun.tif"         "cwd2010jun.tif.aux.xml" "cwd2010jun.tif.xml"    
## [22] "cwd2010mar.tif"         "cwd2010may.tif"         "cwd2010nov.tif"        
## [25] "cwd2010oct.tif"         "cwd2010sep.tif"

almost there… we still have problems with .tif.* files like tif.xml or tif.aux.xml. We can limit these using what I will call an ANTI-wildcard. “$”. '.tif\$' will tell the computer the file needs to END in .tif. Let’s also store the path/file name by setting full.names=TRUE

list.files(path="./data/raster/CWD/",pattern = '.tif$', full.names=TRUE)

##  [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
##  [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
##  [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
##  [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
##  [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"

Let’s store that information:

CWDS = list.files(path="./data/raster/CWD/", pattern = 'tif$', full.names=TRUE)

Challenge: In teams read in and sum all of these rasters, then divide by the total number of files.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

Finally, let’s write out our results:

# writeRaster(x=yourraster,filename="./output/CWD_mn_2010.tif")

Raster Stacks

R also has the capactiy to handle multi-band or multi-layer raster ‘stacks’. In this case we will be using a stack to combine multiple files of the same type (CWD) over time (months of 2010). As you will see there are multiple advantages to handling data in stacks.

The only requirement is that all the images must have the same: - Image extent - Resolution - Projection

In this case our CWD data is already identical in shape and projection. All we need to do is use a list of files like CWDS that we created earlier.

CWDS

##  [1] "./data/raster/CWD//cwd2010apr.tif" "./data/raster/CWD//cwd2010aug.tif"
##  [3] "./data/raster/CWD//cwd2010dec.tif" "./data/raster/CWD//cwd2010feb.tif"
##  [5] "./data/raster/CWD//cwd2010jan.tif" "./data/raster/CWD//cwd2010jul.tif"
##  [7] "./data/raster/CWD//cwd2010jun.tif" "./data/raster/CWD//cwd2010mar.tif"
##  [9] "./data/raster/CWD//cwd2010may.tif" "./data/raster/CWD//cwd2010nov.tif"
## [11] "./data/raster/CWD//cwd2010oct.tif" "./data/raster/CWD//cwd2010sep.tif"

and create a stack from the data:

cwd_stack = stack(CWDS)

Now looking at the properties we can see that cwd_stack has 12 layers.

cwd_stack

## class      : RasterStack 
## dimensions : 1120, 872, 976640, 12  (nrow, ncol, ncell, nlayers)
## resolution : 1080, 1080  (x, y)
## extent     : -375034, 566726, -616284.9, 593315.1  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=0 +lon_0=-120 +lat_1=34 +lat_2=40.5 +x_0=0 +y_0=-4000000 +datum=NAD83 +units=m +no_defs 
## names      : cwd2010apr, cwd2010aug, cwd2010dec, cwd2010feb, cwd2010jan, cwd2010jul, cwd2010jun, cwd2010mar, cwd2010may, cwd2010nov, cwd2010oct, cwd2010sep 
## min values :          0,          0,          0,          0,          0,          0,          0,          ?,          ?,          ?,          ?,          ? 
## max values :  150.96939,  205.35623,   47.32375,   65.93313,   32.79562,  233.96126,  219.87500,          ?,          ?,          ?,          ?,          ?

This allows us to do operations much easier. For instance we can just run arithmetic functions on the entire stack:

CWD_mn_2010 = mean(cwd_stack)
plot(CWD_mn_2010)

We can also do ‘global’ statistics like finding the ‘mean’ of all cells for any given month.

cellStats(cwd_stack, stat='mean', na.rm=TRUE)

##  cwd2010apr  cwd2010aug  cwd2010dec  cwd2010feb  cwd2010jan  cwd2010jul 
##  45.5497388 157.9214946   1.6478612   5.2903873   0.1809682 176.2727036 
##  cwd2010jun  cwd2010mar  cwd2010may  cwd2010nov  cwd2010oct  cwd2010sep 
## 143.6983408  27.6639933  90.9696350  19.2020629  27.7737959 120.8165948

Or other summary visualization like histograms on a month by month basis.

hist(cwd_stack)

Raster Vector Interactions

Raster cropping

Many geographic data projects involve integrating data from many different sources, such as remote sensing images (rasters) and administrative boundaries (vectors). Often the extent of input raster datasets is larger than the area of interest. In this case raster cropping and masking are useful for unifying the spatial extent of input data. Both operations reduce object memory use and associated computational resources for subsequent analysis steps, and may be a necessary preprocessing step before creating attractive maps involving raster data.

We will use two objects to illustrate raster cropping:

• A SpatRaster object srtm representing elevation (meters above sea level) in south-western Utah
• A vector (sf) object zion representing Zion National Park

Both target and cropping objects must have the same projection. The following code chunk therefore not only reads the datasets from the spDataLarge package. This package needs to be installed directly from github using the following code:

install.packages("spDataLarge", repos = "https://geocompr.r-universe.dev")

We will also need to import the sf package to handle vector data.

if(!require("spDataLarge")) install.packages("spDataLarge")

## Loading required package: spDataLarge

if(!require("sf")) install.packages("sf")

srtm = raster(system.file("raster/srtm.tif", package = "spDataLarge"))
zion = read_sf(system.file("vector/zion.gpkg", package = "spDataLarge"))
zion = st_transform(zion, crs(srtm))

Let’s take a look at the two files, note the use of add=TRUE in the second plot command. This adds the new object to the previous plot.

plot(srtm)
plot(zion, add=TRUE)

We use crop() from the raster package to crop the srtm raster. The function reduces the rectangular extent of the object passed to its first argument based on the extent of the object passed to its second argument.

srtm_cropped = crop(srtm, zion)
plot(srtm_cropped)
plot(zion, add=TRUE)

Note that the extent of the strm raster has been cropped to fit the zion polygon.

Related to crop() is the terra function mask(), which sets values outside of the bounds of the object passed to its second argument to NA. The following command therefore masks every cell outside of the Zion National Park boundaries:

srtm_masked = mask(srtm, zion)
plot(srtm_masked)

Importantly, we want to use both crop() and mask() together in most cases. This combination of functions would (a) limit the raster’s extent to our area of interest and then (b) replace all of the values outside of the area to NA.

srtm_cropped = crop(srtm, zion)
srtm_final = mask(srtm_cropped, zion)
plot(srtm_final)

Changing the settings of mask() yields different results. Setting updatevalue = 0, for example, will set all pixels outside the national park to 0. Setting inverse = TRUE will mask everything inside the bounds of the park (see ?mask for details)

srtm_inv_masked = mask(srtm, zion, inverse = TRUE)
plot(srtm_inv_masked)

Raster extraction

Raster extraction is the process of identifying and returning the values associated with a ‘target’ raster at specific locations, based on a (typically vector) geographic ‘selector’ object. The results depend on the type of selector used (points, lines or polygons) and arguments passed to the extract() function, which we use to demonstrate raster extraction. The reverse of raster extraction — assigning raster cell values based on vector objects — is rasterization.

The basic example is of extracting the value of a raster cell at specific points. For this purpose, we will use zion_points, which contain a sample of 30 locations within the Zion National Park.

data("zion_points", package = "spDataLarge")
plot(srtm)
plot(zion, add=TRUE, col=NA)
plot(zion_points, add=TRUE, col='red')

The following command extracts elevation values from srtm and creates a data frame with points IDs (one value per vector’s row) and related srtm values for each point. Now, we can add the resulting object to our zion_points dataset with the cbind() function:

elevation = extract(srtm, zion_points)
zion_points = cbind(zion_points, elevation) # insert the elevation data into the orginal shapefile
zion_points

## Simple feature collection with 30 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2077 ymin: 37.16632 xmax: -112.8717 ymax: 37.43165
## Geodetic CRS:  WGS 84
## First 10 features:
##    elevation                   geometry
## 1       1802 POINT (-112.9159 37.20013)
## 2       2433 POINT (-113.0937 37.39263)
## 3       1886 POINT (-113.0246 37.33466)
## 4       1370 POINT (-112.9611 37.24326)
## 5       1452 POINT (-112.9898 37.20847)
## 6       1635 POINT (-112.8807 37.19319)
## 7       1380 POINT (-113.0505 37.24061)
## 8       2032 POINT (-113.0953 37.34965)
## 9       1830 POINT (-113.0362 37.31429)
## 10      1860 POINT (-113.2077 37.43165)

Note that now we have the elevation values for each one of the points associated with zion_points. For fun let’s use mapview to make an interactive map of both. Again we will use zcol to set the column used for symbology in the point map.

if(!require("mapview")) install.packages("mapview")

## Loading required package: mapview

mapview(srtm)+
mapview(zion_points, zcol='elevation')

Let’s tie this lesson back to our work with OSM. Here we are going to import the package osmdata to allow us to interact with the Overpass Turbo API. Let’s try to build a query using a bounding box from getbb(), where we search for Zion National Park, this bounding box then get passed to Overpass Turbo and restaurants will be added.

if(!require("osmdata")) install.packages("osmdata")

query = getbb(place_name="Zion National Park, Utah, USA") %>%
      opq() %>%
       add_osm_feature(key="amenity", value="restaurant")

restaurant = osmdata_sf(q=query)

Ok, bad request! Thats not too specific. Let’s try to see what happened. First let’s see if we were able to get a bounding box for Zion National Park.

getbb(place_name="Zion National Park, Utah, USA")

Looks like that is the problem. You can try to see if you can adjust the name to get a valid response. But its likely we will need to create a bounding box manually. Luckily we can get the bounding box from the srtm raster file using bbox(srtm), this then needs to be converted to a vector using c() in order to meet the requirement of the opq function. Look at ?opq, notice that we can pass a bounding box in the format bbox = c(xmin, ymin, xmax, ymax).

Note: bounding boxes can also be found using tools like bboxfinder.com.

Now let’s try again, using our new bounding box.

query =  opq(bbox = c(bbox(srtm))) %>%
       add_osm_feature(key="amenity", value="restaurant")

restaurant = osmdata_sf(q=query)
restaurant

## Object of class 'osmdata' with:
##                  $bbox : 37.1320834298579,-113.239583212784,37.5129167631658,-112.85208321281
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 32 points
##             $osm_lines : NULL
##          $osm_polygons : 'sf' Simple Features Collection with 4 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

Now let’s plot it out and assign the restaurant names as labels:

mapview(srtm)+
mapview(restaurant$osm_points,
        label = restaurant$osm_points$name,
        legend=F)

Raster extraction also works with line selectors. Then, it extracts one value for each raster cell touched by a line. In this case, the best approach is to split the line into many points and then extract the values for these points. To demonstrate this, the code below creates zion_transect, a straight line going from northwest to southeast of the Zion National Park

Let’s create a line based on two points c(-113.2, -112.9) and c(37.45, 37.2) which are xy pairs in lat lon. Note that we are assigning a projection using an EPSG code of 4326. EPSG codes are unique idenitifying codes that are availabel for all common projections. To look up EPSG codes please go to https://epsg.io.

zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) %>%
  st_linestring() %>%
  st_sfc(crs =  'EPSG:4326')%>%
  st_sf()
zion_transect

## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
##                         geometry
## 1 LINESTRING (-113.2 37.45, -...

Now let’s see if its showing up in the right place

mapview(srtm)+
mapview(zion_transect, color='white')

The utility of extracting heights from a linear selector is illustrated by imagining that you are planning a hike. The method demonstrated below provides an ‘elevation profile’ of the route (the line does not need to be straight), useful for estimating how long it will take due to long climbs.

Because our line is made up of only two points and extract works best with points, we need to break our longer line into lots of points.

The first step is to add a unique id for each transect. Next, with the st_segmentize() function we can add points along our line(s) with a provided density (dfMaxLength) and convert them into points with st_cast().

Let’s look at st_segmentize help to see what units dfMaxLength is in:

?st_segmentize

Now lets create our collection of points:

zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")
zion_transect

## Simple feature collection with 257 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
## First 10 features:
##     id                   geometry
## 1    1       POINT (-113.2 37.45)
## 1.1  1 POINT (-113.1988 37.44902)
## 1.2  1 POINT (-113.1976 37.44805)
## 1.3  1 POINT (-113.1965 37.44707)
## 1.4  1  POINT (-113.1953 37.4461)
## 1.5  1 POINT (-113.1941 37.44512)
## 1.6  1 POINT (-113.1929 37.44415)
## 1.7  1 POINT (-113.1918 37.44317)
## 1.8  1  POINT (-113.1906 37.4422)
## 1.9  1 POINT (-113.1894 37.44122)

Now, we have a large set of points, and we want to derive a distance between the first point in our transects and each of the subsequent points. In this case, we only have one transect, but the code, in principle, should work on any number of transects:

zion_transect = zion_transect %>%
  group_by(id) %>%
  mutate(dist = st_distance(geometry)[, 1]) 
zion_transect

## Simple feature collection with 257 features and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
## # A tibble: 257 × 3
## # Groups:   id [1]
##       id             geometry  dist
##  * <int>          <POINT [°]>   [m]
##  1     1       (-113.2 37.45)    0 
##  2     1 (-113.1988 37.44902)  150.
##  3     1 (-113.1976 37.44805)  300.
##  4     1 (-113.1965 37.44707)  450.
##  5     1  (-113.1953 37.4461)  600.
##  6     1 (-113.1941 37.44512)  750.
##  7     1 (-113.1929 37.44415)  901.
##  8     1 (-113.1918 37.44317) 1051.
##  9     1  (-113.1906 37.4422) 1201.
## 10     1 (-113.1894 37.44122) 1351.
## # … with 247 more rows

Finally, we can extract elevation values for each point in our transects and combine this information with our main object.

zion_elev = extract(srtm, zion_transect)
zion_transect = cbind(zion_transect, zion_elev)
zion_transect

## Simple feature collection with 257 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -113.2 ymin: 37.2 xmax: -112.9 ymax: 37.45
## Geodetic CRS:  WGS 84
## First 10 features:
##    id          dist zion_elev                   geometry
## 1   1    0.0000 [m]      2001       POINT (-113.2 37.45)
## 2   1  150.0961 [m]      2033 POINT (-113.1988 37.44902)
## 3   1  300.1922 [m]      2020 POINT (-113.1976 37.44805)
## 4   1  450.2883 [m]      1989 POINT (-113.1965 37.44707)
## 5   1  600.3844 [m]      1902  POINT (-113.1953 37.4461)
## 6   1  750.4805 [m]      1869 POINT (-113.1941 37.44512)
## 7   1  900.5766 [m]      1827 POINT (-113.1929 37.44415)
## 8   1 1050.6728 [m]      1767 POINT (-113.1918 37.44317)
## 9   1 1200.7689 [m]      1772  POINT (-113.1906 37.4422)
## 10  1 1350.8650 [m]      1722 POINT (-113.1894 37.44122)

We can now see that we have both the distance from the first point in our line, and the elevation of each point along that line.

Let’s plot the transect and see what it looks like:

plot(x=zion_transect$dist, 
     y = zion_transect$zion_elev, 
     type='l', 
     col='darkred',
     main='Elevation Profile',
     xlab = 'Distance',
     ylab = 'Elevation (m)')

Step 1 - Create a query to pull all the ‘primary’ roads in Zion National Park and save them in a variable called roads. Make sure to look at the tagging convention on OSM’s map features directory. Also you need to set the bbox from the query, you will want to do it based on the srtm raster as we did earlier.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# # hint
# query =  opq(bbox =  ?????) %>%
#        add_osm_feature(key="?????", value="?????")
# 
# roads = osmdata_sf(q=query)

Step 2: osmdata_sf returns points (roads$osm_points), lines (roads$osm_lines), and polygons(roads$osm_polygons). Create a new variable called primary_roads that holds the road line features.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# print(primary_roads)

Step 3: Use the filter() function from dpylr to isolate only roads with the name Zion Park Boulevard, and store the output in a variable called park_road:

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# print(park_road)

Let’s take a look at the data:

park_road

Notice that there are multiple segments making up this one road. Roads are typically broken into segments to allow for connectivity between other roads, and to maintain data like speed limits - which likely vary along the road.

Step 4: We need to ‘dissolve’ all of these road features into a single line representing ‘Zion Park Boulevard’. In order to do this we are going to tell the computer to create a new geometry column using summarize, based on the st_union of the multiple geometry rows. This is one that we haven’t talked about yet, so I will just do it for you.

park_road = park_road %>% 
            summarize(geometry = st_union(geometry))
plot(park_road)

print(park_road)

## Simple feature collection with 1 feature and 0 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -113.0368 ymin: 37.16021 xmax: -112.9902 ymax: 37.19982
## Geodetic CRS:  WGS 84
##                         geometry
## 1 LINESTRING (-113.0368 37.16...

Notice there is now only one row of data and our multiple road segments are now the union.

Step 5: Now we can create our column called id, then use st_segmentize to break the road into multiple points at a fixed distance (250 or 500 works well), then we need to cast that point object back into an sf point object using st_cast.

See if you can find the correct chunk of code

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# park_road$id = 1:nrow(???)
# park_road = st_segmentize(park_road, dfMaxLength = ????)
# park_road = st_cast(park_road, "?????")
# park_road

Step 6: Let’s calculate the distances from the first point to all following points:

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# park_road = park_road %>%
#   group_by(?????) %>%
#   mutate(dist = st_distance(?????)[, 1]) 
# park_road

Step 7: Now we can use the extract function to grab the values from the raster srtm at each one of the points making up park_road, then we need to cbind those values back into the park_road_elev dataset.

### ### ### ### ### 
# YOUR CODE HERE
### ### ### ### ### 

# park_road_elev = extract(????, park_road)
# park_road_elev = cbind(????, ????)
# park_road_elev

Step 8: Plot the park_road_elev$dist on the x axis and park_road_elev$park_road_elev on the y axis, and set the type to 'line'.